% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Warranty Disclaimer and Copyright Notice
% 
% Copyright (C) 2003-2010 Institute for Systems Biology, Seattle, Washington, USA.
% 
% The Institute for Systems Biology and the authors make no representation about the suitability or accuracy of this software for any purpose, and makes no warranties, either express or implied, including merchantability and fitness for a particular purpose or that the use of this software will not infringe any third party patents, copyrights, trademarks, or other rights. The software is provided "as is". The Institute for Systems Biology and the authors disclaim any liability stemming from the use of this software. This software is provided to enhance knowledge and encourage progress in the scientific community. 
% 
% This is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version.
% 
% You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [REGR] = doregress_constrained(DATA,Loc,Weight,val1,val2,opt,val3,constrained_flag,rm);


%% Options
lo = 1;                                                            %log offset (default = 1)
plot_flag = 0;                                                     %plot results (default = 0=>no plotting; 1=>plotting)
if nargin<9;
    % rm is the design matrix for the regression model
    rm(1,:) = [1  1 1 2  2 2 3  3];                                    %regression model (1=FSC, 2=SSC,3=FSC*SSC (interaction effect))
    rm(2,:) = [1 .5 2 1 .5 2 1 .5];                                    %regression model (indicates power to which each predictor is raised)
end
if nargin<8;
    % constrained flag==1=> monotonicity constraints   constrained flag==0=> no monotonicity constraints
    constrained_flag = 1;
end
if nargin<7;
    % when val3 is not empty, but contains two values, additional constraints are added to make sure that the regression surface is between these values.
    % the val3 output of FC_preprocess is an appropriate choice
    val3 = [];
end
if nargin<6;
    % transformation of the data (as explained below)
    opt = 'lin';
end

switch opt %%Transformation of Loc
    case 'lin'  %% Regression without data transformation
        %do nothing
    case 'log'  %% Regression with log data transformation for all three channels (ss,fs,fu)
        Loc = log(Loc+lo);
        val1 = log(val1+lo);
        val2 = log(val2+lo);
    case 'logF'  %% Regression with log data transformation for only the fluoresence channel
        %do nothing
end

N = length(DATA);

for n = 1:N;
    display([num2str(n) ' of ' num2str(N)]);
    D = DATA{n};
    switch opt %%Transformation of data
        case 'lin'  %% Regression without data transformation
            %do nothing
        case 'log'  %% Regression with log data transformation for all three channels (ss,fs,fu)
            D = log(D+lo);
        case 'logF'  %% Regression with log data transformation for only the fluoresence channel
            D(:,3) = log(D(:,3)+lo);
    end
    
    %make regression data;
    Y = D(:,3);
    D(:,3) = D(:,1).*D(:,2);
    X = zeros(length(Y),size(rm,2));
    for m = 1:size(rm,2);
        X(:,m) = D(:,rm(1,m)).^rm(2,m);
    end
    
    B = regress(Y,[ones(size(X,1),1) X]);
    OFFSET = B(1);
    B = B(2:end);
    residual = Y-(X*B+OFFSET);
    
    if plot_flag
        close all; h = figure;hold on;
        plot3(D(:,1),D(:,2),Y,'k.','MarkerSize',4);hold on;
        S = 100;
        XXp = zeros(S,S,5);
        [XXp(:,:,1),XXp(:,:,2)] = meshgrid(linspace(val1(1),val1(2),S),linspace(val2(1),val2(2),S));
        XXp(:,:,3) = XXp(:,:,1).*XXp(:,:,2);
        XXp(:,:,4) = OFFSET;
        for m = 1:size(rm,2);
            XXp(:,:,4) = XXp(:,:,4)+ B(m).*XXp(:,:,rm(1,m)).^rm(2,m);
        end
        hm = mesh(XXp(:,:,1),XXp(:,:,2),XXp(:,:,4));
        hs = surf(XXp(:,:,1),XXp(:,:,2),XXp(:,:,4));
        set(hs,'EdgeColor','none','FaceAlpha',.5);
        set(gca,'Zlim',[0 1.5*max(max(XXp(:,:,4)))]);
        grid;
        set(h,'Position',[1281 1 1280 957]);
        view(-20,28)
        pause;
        set(hs,'FaceAlpha',.1);
        set(hm,'EdgeColor',[.7 .7 .7])
    end
    
    %% Making constraints
    
    if constrained_flag
        
        %1. First choose direction of monotonicity
        S = 5;
        area = .5;
        [bandwidth,P,XX,YY]=kde2d(D(:,[1 2]),2^S,[floor(val1(1)) floor(val2(1))],[ceil(val1(2)) ceil(val2(2))]);
        P = P./sum(sum(P));
        %find 95% of density (or whatever value is in area)
        t = fminsearch(@(t) abs(sum(P(P>t))-area),1/(2^S)^2);
        % imagesc(P>t)
        ZZ = zeros(sum(sum(P>t)),4);
        ZZ(:,1) = XX(P>t);
        ZZ(:,2) = YY(P>t);
        ZZ(:,3) = ZZ(:,1).*ZZ(:,2);
        ZZ(:,4) = OFFSET;
        for m = 1:size(rm,2);
            ZZ(:,4) = ZZ(:,4)+ B(m).*ZZ(:,rm(1,m)).^rm(2,m);
        end
        dir1 = regress(ZZ(:,4),[ones(size(ZZ(:,1))) ZZ(:,1)]);
        dir2 = regress(ZZ(:,4),[ones(size(ZZ(:,1))) ZZ(:,2)]);
        MonoDir = double([dir1(2)>=0 dir2(2)>=0]);
        
        %2. Select points to evaluate on monotonicity
        S = 3;
        [bandwidth,P,XX,YY]=kde2d(D(:,[1 2]),2^S,[floor(val1(1)) floor(val2(1))],[ceil(val1(2)) ceil(val2(2))]);
        noc = 2*(2^S-1)*2^S;
        Q = zeros(noc,4);
        p = 0;
        for i = 1:2^S
            for j = 1:2^S-1
                p = p + 1;
                if MonoDir(1)==1;
                    Q(p,:) = [XX(1,j) YY(i,1) XX(1,j+1) YY(i,1)];
                else
                    Q(p,:) = [XX(1,j+1) YY(i,1) XX(1,j) YY(i,1)];
                end
            end
        end
        for i = 1:2^S
            for j = 1:2^S-1
                p = p + 1;
                if MonoDir(2)==1;
                    Q(p,:) = [XX(1,i) YY(j,1) XX(1,i) YY(j+1,1)];
                else
                    Q(p,:) = [XX(1,i) YY(j+1,1) XX(1,i) YY(j,1)];
                end
            end
        end
        XX1 = zeros(noc,3);
        XX1(:,1) = Q(:,1);
        XX1(:,2) = Q(:,2);
        XX1(:,3) = XX1(:,1).*XX1(:,2);
        A1 = zeros(noc,size(rm,2));
        for m = 1:size(rm,2);
            A1(:,m) = XX1(:,rm(1,m)).^rm(2,m);
        end
        XX2 = zeros(noc,3);
        XX2(:,1) = Q(:,3);
        XX2(:,2) = Q(:,4);
        XX2(:,3) = XX2(:,1).*XX2(:,2);
        A2 = zeros(noc,size(rm,2));
        for m = 1:size(rm,2);
            A2(:,m) = XX2(:,rm(1,m)).^rm(2,m);
        end
        
        if ~isempty(val3);
            
            A = zeros(noc+2,size(rm,2));
            A(1:noc,:) = A1-A2;
            
            XBP = zeros(2,3);
            if MonoDir(1)==1&MonoDir(2)==1;
                XBP(:,1) = XX(1,[1 2^S]);
                XBP(:,2) = YY([1 2^S],1);
            elseif MonoDir(1)==0&MonoDir(2)==1;
                XBP(:,1) = XX(1,[2^S 1]);
                XBP(:,2) = YY([1 2^S],1);
            elseif MonoDir(1)==1&MonoDir(2)==0;
                XBP(:,1) = XX(1,[1 2^S]);
                XBP(:,2) = YY([2^S 1],1);
            elseif MonoDir(1)==0&MonoDir(2)==0;
                XBP(:,1) = XX(1,[2^S 1]);
                XBP(:,2) = YY([2^S 1],1);
            end
            XBP(:,3) = XBP(:,1).*XBP(:,2);
            AX = zeros(2,size(rm,2));
            for m = 1:size(rm,2);
                AX(:,m) = XBP(:,rm(1,m)).^rm(2,m);
            end
            
            A(noc+1,:) = -AX(1,:);
            A(noc+2,:) = AX(2,:);
            
            b = zeros(noc+2,1);
            b(noc+1) = val3(1);
            b(noc+2) = val3(2);
            noc = noc + 2;
            
        else
            
            A = A1-A2;
            b = zeros(noc,1);
            
        end
        
        %constrained regression analysis
        options = optimset('LargeScale','off');
        NB = zeros(size(B));
        [B2,resnorm,Nresidual,exitflag,output,lambda] = lsqlin([X ones(size(Y))],Y,[A zeros(noc,1)],b,[],[],[],[],[],options);
        NOFFSET = B2(end);
        NB = B2(1:end-1);
        Nresidual = -Nresidual;
        
    else
        NB = B;
        NOFFSET = OFFSET;
        Nresidual = residual;
    end
    
    if plot_flag
        XXp(:,:,5) = NOFFSET;
        for m = 1:size(rm,2);
            XXp(:,:,5) = XXp(:,:,5)+ NB(m).*XXp(:,:,rm(1,m)).^rm(2,m);
        end
        hm = mesh(XXp(:,:,1),XXp(:,:,2),XXp(:,:,5));
        hs = surf(XXp(:,:,1),XXp(:,:,2),XXp(:,:,5));
        set(hs,'EdgeColor','none','FaceAlpha',.5);
        pause
    end
    
    %Offset area
    XO = zeros(length(Weight),4);
    XO(:,1) = Loc(:,1);
    XO(:,2) = Loc(:,2);
    XO(:,3) = XO(:,1).*XO(:,2);
    XO(:,4) = NOFFSET;
    for m = 1:size(rm,2);
        XO(:,4) = XO(:,4)+ NB(m).*XO(:,rm(1,m)).^rm(2,m);
    end
    
    R1 = Nresidual + ((Weight'*XO(:,4))/sum(Weight));
    
    switch opt %%Results for REGR
        case 'lin'  %% Regression without data transformation
            REGR{n} = R1;
        case 'log'  %% Regression with log data transformation for all three channels (ss,fs,fu)
            REGR{n} = exp(R1)-lo;
        case 'logF'  %% Regression with log data transformation for only the fluoresence channel
            REGR{n} = exp(R1)-lo;
    end
    
end




